---
  title: "Winbugs for Bayesian mediation analysis"
output:
  pdf_document: default
word_document: default
html_document: default
---
  
  ```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

To use the R package R2WinBUGS, we first install the package in R (`install.packages("R2WinBUGS")`) and load it. 
```{r}
library(R2WinBUGS)
```

#When the mediator is binary
##Method 1
the usual method 1 does not work since we cannot use alpha directly.

##Method 2
generate the data:
```{r}
mu_m=alpha*x
M3=rbinom(100,1,exp(mu_m)/(1+exp(mu_m)))
y2.1=c*x+beta*M3+e2

```

##Method 2
The product of partial difference still wookrs. The only change is in the model: need to set up a link function and random term for the third variable.
```{r}
deltax=diff(range(x))/10
deltam=diff(range(M))/10
data1<- list ( x=x,M=M3,y=y2.1,N=N,deltax=deltax,deltam=deltam)
inits<- function(){list(alpha=0,c=0,beta=0,prec2=0.5)}
med3.2<- bugs(data1, inits, model.file = "O:/My Documents/My Research/Research/Bayesian Mediation Analysis/Winbugs/mediation_c_b_c_2.txt",
               parameters = c("alpha","c","beta","var2",
                                 "ie","de","te"),
               n.chains = 1, n.iter = 11000,n.burnin=1000,
               bugs.directory = "O:/Current work/Program/winbugs143_unrestricted/winbugs14_full_patched/WinBUGS14/",
               debug = F)

mean(med3.2$mean$alpha)
mean(med3.2$mean$beta)
mean(med3.2$mean$c)
mean(med3.2$mean$ie)
mean(med3.2$mean$de)
mean(med3.2$mean$te)
```

##Method 3
Use the marginal distribution of M
```{r}
data1<- list ( x=x,M=M3,y=y2.1,N=N,deltax=deltax)
inits<- function(){list(alpha=0,c=0,beta=0,prec1=0.5)}
med3.3<- bugs(data1, inits, model.file = "O:/My Documents/My Research/Research/Bayesian Mediation Analysis/Winbugs/mediation_c_b_c_3.txt",
               parameters = c("alpha","c","beta","var2",
                       "ie","de","te"),
               n.chains = 1, n.iter = 15000,n.burnin=5000,
               bugs.directory = "O:/Current work/Program/winbugs143_unrestricted/winbugs14_full_patched/WinBUGS14/",
               debug = F)

mean(med3.3$mean$alpha)
mean(med3.3$mean$beta)
mean(med3.3$mean$c)
mean(med3.3$mean$ie)
mean(med3.3$mean$de)
mean(med3.3$mean$te)

```

In summary, when the mediator is binary, the main change is to set up m: need to put a link on (mu_M) and use a different random component.

#When the mediator is binary
In the model setting, we set
    log(mu_M[i]) <-  c*x[i]+beta*M[i]
    y[i] ~ dpois(mu_y[i])
In general, for different outcome, we change the link function and the random term.